home *** CD-ROM | disk | FTP | other *** search
/ PC Electronics Plus 3 / PC Electronics Plus 3.iso / coil01 / formulae.c < prev    next >
C/C++ Source or Header  |  1993-02-11  |  15KB  |  730 lines
  1.  
  2. /* NBS inductance formulas
  3.  *
  4.  *  Program by Steve Moshier
  5.  *  moshier@world.std.com
  6.  *
  7.  * Last rev: December, 1992
  8.  */
  9.  
  10.  
  11. #define DEBUG 0
  12.  
  13. #define SQRT2 1.41421356237309504880
  14. #define PI 3.14159265358979323846
  15. #define fabs(x) ( (x) < 0 ? -(x) : (x) )
  16.  
  17. double MACHEP = 2.2e-16;
  18. double MAXNUM = 1.7e38;
  19. double Kn, Kp, P, kl;
  20. double sqrt(), log(), pow(), sin(), cos(), atan();
  21. double ellpe(), ellpk(), chbevl();
  22. double rsol(), lyle(), dwighta(), dwightb(), butterw();
  23. double lyle2(), spielrein(), nagaoka();
  24.  
  25.  
  26. /* Dispatch program for circular coil
  27.  * of arbitrary rectangular winding cross section
  28.  * decides which formula to use.
  29.  */
  30. double rsol( a, b, c, N )
  31. double a, b, c, N;
  32. {
  33. double r, s, t, gam, c2a, bc, b2a, L;
  34.  
  35. if( c == 0.0 )
  36.     {
  37.     L = nagaoka( a, b, N );
  38. #if DEBUG
  39.     printf( "N" );
  40. #endif
  41.     goto soldone;
  42.     }
  43. bc = b/c;
  44. c2a = 0.5*c/a;
  45. b2a = 0.5*b/a;
  46. r = a + 0.5*c;
  47. s = a - 0.5*c;
  48. gam = s/r;
  49.  
  50. if( b == 0.0 )
  51.     {
  52.     if( gam < 0.5 )
  53.         {
  54. #if DEBUG
  55.         printf( "SP" );
  56. #endif
  57.         L = spielrein( a, b, c, N );
  58.         goto soldone;
  59.         }
  60.     else
  61.         {
  62. #if DEBUG
  63.         printf( "L2P" );
  64. #endif
  65.         L = lyle2( a, b, c, N );
  66.         goto soldone;
  67.         }
  68.     }
  69.  
  70. t = c2a*c2a * (1.0 + bc*bc);
  71. if( t <= 1.0 )
  72.     {
  73. #if DEBUG
  74.     printf( "LP" );
  75. #endif
  76.     L = lyle( a, b, c, N );
  77.     goto soldone;
  78.     }
  79. if( (gam <= 0.4) && (r/b <= 0.4) )
  80.     {
  81. #if DEBUG
  82.     printf( "DKp" );
  83. #endif
  84.     L = dwightb( a, b, c, N );
  85.     goto soldone;
  86.     }
  87. if( (c2a <= 0.3) && (b2a > 0.7) )
  88.     {
  89. #if DEBUG
  90.     printf( "Dk" );
  91. #endif
  92.     L = dwighta( a, b, c, N );
  93.     goto soldone;
  94.     }
  95.  
  96. #if DEBUG
  97. printf( "BKp" );
  98. #endif
  99. L = butterw( a, b, c, N );
  100. soldone:
  101. return(L);
  102. }
  103.  
  104. double log(), sqrt(), asin(), atan();
  105. extern double fac[]; /* factorial table */
  106.  
  107.  
  108. /* Circular solenoidal current sheet (infinitely thin tape winding)
  109.  *
  110.  * Nagaoka (1909): J. College of Science, Tokyo, 27, Art. 6, p. 18
  111.  *
  112.  * See also E. Jahnke and F. Emde, _Tables of Functions_,
  113.  * 4th edition, Dover, 1945, pp 86-89.
  114.  */
  115. double nagaoka( a, b, N )
  116. double a, b, N;
  117. {
  118. double t, f, k1e, ke, Ke, Ee, tana, Ls;
  119.  
  120.  
  121. if( a <= 0.0 )
  122.     {
  123.     goto tiny;
  124.     }
  125. if( b <= 0.0 )
  126.     {
  127.     if( a <= 0.0 )
  128.         {
  129.         Kn = 0.0;
  130.         Kp = 0.0;
  131.         return( 0.0 );
  132.         }
  133.     else
  134.         {
  135. tiny:
  136.         Kn = 1.0;
  137.         Kp = 1.0;
  138.         return( MAXNUM );
  139.         }
  140.     }
  141.  
  142. tana = 2.0*a/b;        /* form factor */
  143. if( tana < 1.0e-7 )
  144.     {
  145.     Kn = 1.0;
  146.     Kp = 1.0;
  147.     P = 4.0*PI*PI*a*Kp/b;
  148.     t = 1.0e-9 * N * N * a * P;
  149.     return(t);
  150.     }
  151. t = atan( tana );
  152. ke = sin(t);     /* modulus, k, of the elliptic integrals */
  153. k1e = cos(t);        /* subroutine argument = 1 - k^2 */
  154. t = k1e * k1e;
  155. Ke = ellpk(t);      /* See elliptic integrals below */
  156. Ee = ellpe(t);
  157. t = tana * tana;  /* square of tan(a) */
  158. /* Note that the constant of proportionality is exact.
  159.  * The "permeability of the vacuum" is 4 pi 10^-9 henry/cm.
  160.  */
  161. f = 4.0*PI*((Ke + (t - 1.0) * Ee)/ke  -  t)/3.0;
  162. t = 2.0e-9 * N * N * a * f; /* the inductance */
  163. Kn = f * b / (2.0*PI*PI*a); /* Nagaoka's constant */
  164. Kp = Kn;
  165. P = 4.0*PI*PI*a*Kp/b;
  166. return(t);
  167. }
  168.  
  169.  
  170.  
  171. /* Single-layer circular helical solenoid of round wire
  172.  *
  173.  * c is the diameter of the wire.
  174.  * There are N complete turns evenly spaced over length b.
  175.  * b is measured from wire center at beginning of first
  176.  *   turn to wire center at end of last turn.
  177.  * a is the coil radius measured to the center of the wire.
  178.  *
  179.  * Snow, 1952
  180.  */
  181. double helsol( a, b, c, N )
  182. double a, b, c, N;
  183. {
  184. double t, f, Ke, Ee, ke, k1e, tana, Ls, z;
  185.  
  186. /* first calculate the thin solenoid inductance */
  187. tana = 2.0*a/b;        /* form factor */
  188. t = atan( tana );
  189. ke = sin(t);     /* modulus, k, of the elliptic integrals */
  190. k1e = cos(t);        /* subroutine argument = 1 - k^2 */
  191. t = k1e * k1e;
  192. Ke = ellpk(t);      /* See elliptic integrals below */
  193. Ee = ellpe(t);
  194. t = tana * tana;  /* square of tan(a) */
  195. f = 4.0*PI*((Ke + (t - 1.0) * Ee)/ke  -  t)/3.0;
  196. Ls = 2.0 * N * N * a * f; /* the inductance of the solenoidal sheet */
  197.  
  198. z = PI*N*c/b;
  199. f = 2.0*N*(0.25-log(z))
  200.   + log(2.0*PI*N*a/b)/3.0
  201.   - (4.0/(PI*PI)) * (Ee/ke - 1.0) * (1.0 + z*z/8.0)
  202.   - (2.0/3.0) * ( (Ke-Ee)/ke - 0.5*ke*Ke )
  203.   - (0.5*k1e/ke) * ( 1.0 - (k1e/ke)*asin(ke) );
  204. f *= 2.0*PI*a;
  205. f +=  b * (log((1.0+k1e)/(1.0-k1e)) + k1e*log(4.0));
  206. f += Ls;
  207. return( 1.0e-9*f );
  208. }
  209.  
  210.  
  211. /* Single-layer square solenoid
  212.  * a = length of the side of the square coil form
  213.  * b = length of the winding
  214.  * N = number of turns
  215.  *
  216.  * Grover 1922b
  217.  */
  218. double squaresol( a, b, N )
  219. double a, b, N;
  220. {
  221. double sum, c1, c2, c3, c4, x1, x2;
  222.  
  223. /* short coil */
  224. if( b <= 0.25*a )
  225. {
  226. x1 = b/a;
  227. c1 = (3.0-2.0*SQRT2)/24.0;
  228. c2 = (6.0*SQRT2-5.0)/480.0;
  229. c3 = (23.0*SQRT2-14.0)/5376.0;
  230. c4 = -(37.0*SQRT2-28.0)/9216.0;
  231. sum = -log(x1) + 0.72599 + x1/3.0;
  232. x2 = x1*x1;
  233. sum += (((c4*x2 + c3)*x2 + c2)*x2 + c1)*x2;
  234. goto done;
  235. }
  236.  
  237. /* long coil */
  238. if( a < 0.25*b )
  239. {
  240. x1 = a/b;
  241. x2 = x1*x1;
  242. c1 = 0.473201004409338551642; /* (2/pi)*(log(1+sqrt(2)) - (sqrt(2)-1)/3) */
  243. c2 = 1.0/(2.0*PI);
  244. c3 = -1.0/(12.0*PI);
  245. c4 = 1.0/(28.0*PI);
  246. sum = 1.0 - c1*x1 + ((c4*x2 + c3)*x2 + c2)*x2;
  247. sum *= 0.5*PI*x1;
  248. goto done;
  249. }
  250.  
  251. /* general closed-form solution */
  252. x1 = a/b;
  253. c1 = a*a + b*b;
  254. c2 = sqrt(c1);
  255. c3 = a*a + c1;
  256. c4 = sqrt(c3);
  257. x2 = x1*x1;
  258. sum = log((a+c2)/(a+c4));
  259. sum += x2 * log((a+c4)/((1.0+SQRT2)*a));
  260. sum += 0.5*(1.0-x2) * log(c1);
  261. sum += x2 * log(a);
  262. sum -= log(b);
  263. sum -= PI*x1;
  264. x2 = (2.0*a + c4)/(SQRT2*(a + c4));
  265. sum += 2.0*x1 * asin(x2);
  266. if( x1 == 0.0 )
  267.     sum += PI * x1;
  268. else
  269.     sum += 2.0 * x1 * atan(1/x1);
  270. sum += (c2-c4)*x1/b;
  271. sum += (SQRT2-1.0)*x1*x1/3.0;
  272. sum += (c3*c4 - 2.0*c1*c2)/(3.0*a*b*b);
  273. sum += b/(3.0*a);
  274.  
  275. done:
  276. return( 8.0e-9*a*N*N*sum );
  277. }
  278.  
  279.  
  280. /* Butterworth's expansions for thick circular coils
  281.  * Grover 1922a
  282.  */
  283. double butterw( a, b, c, N )
  284. double a, b, c, N;
  285. {
  286. double r, s, x, y, z, gam, g3, zz, c2a;
  287. double butchi(); /* see below */
  288.  
  289. r = a + 0.5*c; /* outer radius */
  290. s = a - 0.5*c; /* inner radius */
  291. gam = s/r;
  292. z = b/r;
  293. zz = gam * gam;
  294. g3 = gam * zz;
  295. if( (z > 0.0) && (gam < 1.0e-12) )
  296.     x = 0.0;
  297. else
  298.     x =  zz * butchi( z/gam, 1.0 );
  299. x -= 2.0 * butchi( z, gam );
  300. x *= g3;
  301. x += butchi( z, 1.0 );
  302. /* y part */
  303. x -=  (1.0 + zz*g3) * butchi( 0.0, 1.0 );
  304. if( gam <= 0.0 )
  305.     x += g3/3.0;
  306. else
  307.     x += 2.0 * g3 * butchi( 0.0, gam );
  308. c2a = 0.5*c/a;
  309. zz = 1.0 + c2a;
  310. z = 2.0*a/b;
  311. y = (3.*gam - 4.)*g3 + 1.0 + 3.0*z*zz*x;
  312. zz = zz * zz;
  313. Kp = (y*zz*zz)/(24.0*c2a*c2a);
  314. P = 2.0*PI*PI*z*Kp;
  315. /*y = (2.0e-9*PI*PI)*z*N*N*a*Kp;*/
  316. y = 1.0e-9*a*N*N*P;
  317. return(y);
  318. }
  319.  
  320.  
  321. double butchi(z, gam)
  322. double z, gam;
  323. {
  324. double zz, u, v, sum, c1, c2, c3, c4, alpha, beta, m;
  325. double zet, zet2, zetp1, chi;
  326. int n;
  327.  
  328.  
  329. if( z > 4.0 )
  330. {
  331. zz = gam*gam;
  332. c1 = ((((14./33.)*zz + (28./9.))*zz + (40./7.))*zz + (28./9.))*zz + (14./33.);
  333. c2 = (((5./27.)*zz + (6./7.))*zz + (6./7.))*zz + (5./27.);
  334. c3 = ((2./21.)*zz + (6./25.))*zz + (2./21.);
  335. c4 = (1./15.)*zz + (1./15.);
  336. zz = 0.25/(z*z);
  337. chi = (((c1*zz - c2)*zz + c3)*zz - c4)*zz + (1./9.);
  338. chi /= 2.0 * z;
  339. return(chi);
  340. }
  341.  
  342. if( z > 0.0 )
  343. {
  344. zet2 = 1.0 + z * z;
  345. zet = sqrt( zet2 );
  346. zetp1 = 1.0 + zet;
  347. zz = 1.0/zet2;
  348. c1 = ((( -35.*zz + 35.)*zz - 3.)*zz - 1.)*zz + 1.0/(zetp1*zetp1);
  349. c1 *= 3./(82944.*zet);
  350. c2 = ( -3.*zz + 1.)*zz + (1./zetp1);
  351. c2 *=  -1.0/(1344.*zet);
  352. }
  353.  
  354. if( z > gam )
  355. {
  356.  
  357. zz = log(zetp1/z);
  358. c3 = (zz  - 1.0/zet)/40.0;
  359. c4 = (0.5*z*z*zz + 0.5*zet - z)/3.0;
  360.  
  361. if( gam <= 0.0 )
  362.     {
  363.     return( c4 );
  364.     }
  365.  
  366. /*   inf      n
  367.  *    -   (-1)  (2n-3)!            2n
  368.  *    >   ---------------- (gam/2z)
  369.  *    -   (2n+3) n! (n+1)!
  370.  *   n=2
  371.  */
  372. zz = 0.5*gam/z;
  373. zz = zz * zz;
  374. v = zz * zz;
  375. sum = 0.0;
  376. n = 2;
  377. do
  378.     {
  379.     u = v*fac[2*n-3]/((2.*n+3.)*fac[n]*fac[n+1]);
  380.     sum += u;
  381.     v *= -zz;
  382.     n += 1;
  383.     }
  384. while( fabs(u/sum) > 1.0e-10 );
  385.  
  386. alpha = (sum*z*z)/(gam*gam);
  387. zz = gam*gam;
  388. chi = (( -c1*zz + c2)*zz - (c3-alpha))*zz + c4;
  389. return(chi);
  390. }
  391.  
  392. if( z > 0.0 )
  393. {
  394. /*    8         n
  395.  *    -   4 (-1)  (2n-3)!                       2n
  396.  *    >   ------------------------   (z / 2 gam)
  397.  *    -   (n-1)!(n+1)!(2n+1)(2n+3)
  398.  *   n=2
  399.  */
  400. zz = 0.5*z/gam;
  401. zz = zz * zz;
  402. v = zz*zz;
  403. sum = 0.0;
  404. for( n=2; n<=8; n++ )
  405.     {
  406.     sum += v*fac[2*n-3]/(fac[n-1]*fac[n+1]*(2.*n+1.)*(2.*n+3.));
  407.     v *= -zz;
  408.     }
  409. zz = z/gam;
  410. sum = 4.0 * zz * zz * sum;
  411. beta = ((( -(log(2.0/zz) + (77./60.))/30.)*zz + (1./6.))*zz - (1./3.))*zz*zz;
  412. beta = 0.5*zz*( (1./3.) + beta - sum );
  413.  
  414. zz = log( 2.0*zetp1 );
  415. c3 = (zz - (1.0/zet) + (29./20.))/40.;
  416. c4 = ( z*z*(0.5*zz - (1./3.)) - z + 0.5*zet)/3.;
  417. zz = gam*gam;
  418. chi = (( -c1*zz + c2)*zz - (c3-beta))*zz + c4;
  419. chi += (z*z/6. - zz/40.) * log(1./gam);
  420. return(chi);
  421. }
  422.  
  423. /* else z == 0 */
  424.  
  425. if( gam == 1.0 )
  426.     return( 0.12206342136 );
  427. if( gam <= 0.0 )
  428.     return( 1.0/6.0 );
  429.  
  430. /*   inf                  2            2n
  431.  *    -   (1*3*5...(2n+3))          gam
  432.  *    >   (--------------)  ----------------------
  433.  *    -   (2*4*6...(2n+4))  (2n+7)(2n+3)(n+3)(n+1)
  434.  *   n=0
  435.  */
  436. c1 = 0.5;
  437. v = 1.0;
  438. zz = gam*gam;
  439. m = 0.0;
  440. sum = 0.0;
  441. do
  442.     {
  443.     c1 *= (2.*m+3.)/(2.*m+4.);
  444.     u = (c1*c1*v)/((2.*m+7.)*(2.*m+3.)*(m+3.)*(m+1.));
  445.     sum += u;
  446.     v *= zz;
  447.     m += 1.0;
  448.     }
  449. while( fabs(u/sum) > 1.0e-8 );
  450. chi = (1./6.) - zz*(log(4./gam) + (9./20.))/40. + 0.5*zz*zz*sum;
  451. return(chi);
  452. }
  453.  
  454.  
  455. /* Lyle's formumla for narrow disk coils, b=0
  456.  */
  457. double lyle2( a, b, c, N )
  458. double a, b, c, N;
  459. {
  460. double c1, c2, c3, c4, c2a, zz, x;
  461.  
  462. c2a = 0.5*c/a;
  463. x = log( 4.0/c2a );
  464. zz = c2a * c2a;
  465.  
  466. c1 = (103./(105.*1024.)) * (x + 1.1394);
  467. c2 = (11./2880.)*(x + (96./55.));
  468. c3 = (x + (43./12.))/24.0;
  469. c4 = x - 0.5;
  470.  
  471. P = (((c1*zz + c2)*zz + c3)*zz + c4) * (4.0*PI);
  472. Kp = 0.0;
  473. x = 1.0e-9*N*N*a*P;
  474. return(x);
  475. }
  476.  
  477.  
  478. /* Spielrein's formula for wide disk coils, b=0
  479.  *
  480.  * Spielrein (1915): Archiv fur Elektrotechnik 3, p. 182
  481.  */
  482. double spielrein( a, b, c, N )
  483. double a, b, c, N;
  484. {
  485. double gam, zz, S, c2a;
  486.  
  487. gam = (a - 0.5*c)/(a + 0.5*c);
  488. S = 6.96957;
  489. if( gam > 1.0e-12 )
  490.     {
  491.     zz = gam*gam;
  492.     S += ((((0.0337*zz
  493.      + 0.06038)*zz
  494.      + 0.12494)*zz
  495.      + 0.33045)*zz
  496.      + 1.48044)*zz*zz*gam;
  497.  
  498.     S +=  -zz*gam*( 13.1595*log(1./gam) + 9.08008 );
  499.     }
  500. c2a = 0.5*c/a;
  501. zz = 1.0 + c2a;
  502. P = (zz*zz*zz*S)/(4.0*c2a*c2a);
  503. Kp = 0.0;
  504. zz = 1.0e-9*N*N*a*P;
  505. return(zz);
  506. }
  507.  
  508. /* Dwight's formula for long, thin coils
  509.  */
  510. double dwighta( a, b, c, N )
  511. double a, b, c, N;
  512. {
  513. double ab, p, pp, c1, c2, c3, c4, c5;
  514.  
  515. ab = 2.0*a/b;
  516. p = 1.0/sqrt( 1.0 + 4.0/(ab*ab) );
  517. pp = p * p;
  518. c1 = ((( (1183./96.)*pp
  519.  - (1117./672.))*pp
  520.  + (15./112.))*pp
  521.  - (1./120.))*pp*p;
  522.  
  523. c2 = ((((( -(3913./128.)*pp
  524.  + (38857./4608.))*pp
  525.  - (1265/576.))*pp
  526.  + (53./96.))*pp
  527.  - (17./180.))*pp
  528.  + (1./36.))*p;
  529.  
  530. c3 = ((((((( -(13.*68915./32768.)*pp
  531.  + (11.*1961./2048.))*pp
  532.  - (2135./512.))*pp
  533.  + (217./128.))*pp
  534.  - (95./128.))*pp
  535.  + (1./3.))*pp
  536.  - (5./24.))*pp
  537.  + (1./6.))*p;
  538.  
  539. p = 0.5*c/a;
  540. pp = p * p;
  541.  
  542. c4 = (( (4547./(350.*3584.))*pp
  543.  + (1./400.))*pp
  544.  - (23./12.))*pp;
  545.  
  546. c5 = (( -(23./4480.)*pp
  547.  - (1./20.))*pp
  548.  + 1.0)*pp;
  549.  
  550. kl = ((1./(3.*PI))*(c5*log(4./p) + c4)
  551.  + ((c1*pp + c2)*pp + c3)*pp)*ab
  552.  + (1./3.)*pp
  553.  - (2./3.)*p;
  554. kl = -kl;
  555. p = nagaoka( a, b, N );
  556. Kp = Kn - kl;
  557. P = 4.0*PI*PI*Kp*a/b;
  558. p = 1.0e-9*N*N*a*P;
  559. return(p);
  560. }
  561.  
  562.  
  563. /* Dwight's formula for long, thick coils
  564.  * Dwight (1918): Electrical World 71, p. 300
  565.  */
  566. double dwightb( a, b, c, N )
  567. double a, b, c, N;
  568. {
  569. double r, s, gam, u, v, sum, zz;
  570. double g3, g4, g5, g7, g9, Q;
  571. double c1, c2, c3, c4, c5, rb;
  572. int n;
  573.  
  574. r = a + 0.5*c;
  575. s = a - 0.5*c;
  576.  
  577. /* Check for solid coil, inner radius zero */
  578. if( (s < 1.0e-12) && (b > 0.0) )
  579. {
  580. u = 2.0*a/b;
  581. if( u > 0.3 )
  582.     return( butterw( a, b, c, N ) );
  583.  
  584. c1 = 197./2016.;
  585. c2 = 113./1400.;
  586. c3 = 1./10.;
  587. c4 = 1./3.;
  588. c5 = 0.7323816;
  589. v = u * u;
  590. sum = ((( -c1*v + c2)*v - c3)*v + c4)*v - c5 * u + 1.0;
  591. P = 8.0*PI*PI*a*sum/(3.0*b);
  592. Kp = P*b/(4.0*PI*PI*a);
  593. u = 1.0e-9*a*N*N*P;
  594. return(u);
  595. }
  596.  
  597. gam = s/r;
  598. /*   inf
  599.  *    -   (2n-1)! (2n+1)!                  2n
  600.  *    >   ---------------------- (gam / 4 )
  601.  *    -   n!n!(n+1)!(n+2)!(2n+5)
  602.  *   n=2
  603.  */
  604. sum = 0.0;
  605. zz = 0.25*gam;
  606. zz = zz*zz;
  607. v = 1.0;
  608. n = 1;
  609. do
  610.     {
  611.     v *= zz;
  612. u = (v*fac[2*n-1]*fac[2*n+1])/(fac[n]*fac[n]*fac[n+1]*fac[n+2]*(2.*n+5.));
  613.     sum += u;
  614.     n += 1;
  615.     }
  616. while( fabs(u/sum) > 1.0e-7 );
  617.  
  618. zz = gam*gam;
  619. g3 = gam*zz;
  620. g4 = zz*zz;
  621. g5 = g4*gam;
  622. g7 = g4*g3;
  623. g9 = g4*g4*gam;
  624. c1 = (9./112.)*(1.0-g5)*(1.0-g7) + (5./288.)*(1.0-g3)*(1.0-g9);
  625. c2 = (9./200.)*(1.0-g5)*(1.0-g5) + (1./28.)*(1.0-g3)*(1.0-g7);
  626. c3 = (1./10.)*(1.0-g3)*(1.0-g5);
  627. c4 = (1./3.)*(1.0-g3)*(1.0-g3);
  628. rb = r/b;
  629. zz = rb*rb;
  630. Q = ((( -c1*zz + c2)*zz - c3)*zz + c4)*zz + 3.0*g5*sum*rb;
  631. Q = Q - 3.*rb*( 0.2441272 - (2./3.)*g3 + (0.4277559 - 0.1*log(gam))*g5);
  632. Q = Q + 1.0 - 4.0*g3 + 3.0*g3*gam;
  633. u = 0.5*c/a;
  634. v = 1.0 + u;
  635. v = v * v;
  636. Kp = (v*v*Q)/(24.0*u*u);
  637. P = 4.0*PI*PI*Kp*a/b;
  638. u = 1.0e-9*N*N*a*P;
  639. return(u);
  640. }
  641.  
  642.  
  643.  
  644. /* Lyle's formula for short, thin coils
  645.  *
  646.  * Lyle (1914), Phil. Trans. 213A, pp. 421-435
  647.  */
  648. #include "lyle.h"
  649.  
  650. double lyle( a, b, c, N )
  651. double a, b, c, N;
  652. {
  653. double x, d, p, q, m1, m2, m3, l0, l1, l2, l3;
  654. double b2, c2;
  655.  
  656. if( c == 0.0 )
  657.     goto smc;
  658. x = b/c;
  659. if( x <= 1.0 )
  660.     {
  661.     x = 2.0*x - 1.0;
  662. /*
  663.  *    m1 = 1.0e-2*chbevl( x, lm1, 10 );
  664.  *    m2 = 1.0e-4*chbevl( x, lm2, 10 );
  665.  */
  666.     m3 = 1.0e-6*chbevl( x, lm3, 10 );
  667.     l0 = chbevl( x, ll0, 14 );
  668.     l1 = 1.0e-2*chbevl( x, ll1, 12 );
  669.     l2 = 1.0e-4*chbevl( x, ll2, 10 );
  670.     l3 = 1.0e-6*chbevl( x, ll3, 10 );
  671.     }
  672. else
  673.     {
  674. smc:
  675.     x = c/b;
  676.     x = 2.0*x - 1.0;
  677. /*
  678.  *    m1 = 1.0e-2*chbevl( x, lm1a, 10 );
  679.  *    m2 = 1.0e-4*chbevl( x, lm2a, 10 );
  680.  */
  681.     m3 = 1.0e-6*chbevl( x, lm3a, 8 );
  682.     l0 = chbevl( x, ll0a, 14 );
  683.     l1 = 1.0e-2*chbevl( x, ll1a, 10 );
  684.     l2 = 1.0e-4*chbevl( x, ll2a, 10 );
  685.     l3 = 1.0e-6*chbevl( x, ll3a, 10 );
  686.     }
  687. b2 = b * b;
  688. c2 = c * c;
  689. m1 = (1./96.)*(3.0*b2 + c2)/(b2 + c2);
  690. d = b2 + c2;
  691. m2 = (1./92160.)*(-90.*b2*b2 + 105.*b2*c2 + 22.*c2*c2)/(d*d);
  692. d = sqrt(d);
  693. x = d/a;
  694. x = x * x;
  695. p = ((m3*x + m2)*x + m1)*x + 1.0;
  696. q = ((l3*x + l2)*x + l1)*x - l0;
  697. P = 4.0*PI*(p * log( 8.0*a/d )  +  q);
  698. Kp = P*b/(4.0*PI*PI*a);
  699. x = 1.0e-9*a*N*N*P;
  700. return(x);
  701. }
  702.  
  703. /* compute Chebyshev polynomials up to degree n-1
  704.  * at arg x
  705.  * and evaluate Chebyshev expansion.
  706.  * Note, the zero degree term is not multiplied by 1/2.
  707.  */
  708. double chbevl( x, ch, n )
  709. double x;
  710. double ch[];
  711. int n;
  712. {
  713. double t[15];
  714. double s;
  715.  
  716. int i;
  717.  
  718. t[0] = 1.0;
  719. t[1] = x;
  720. s = x + x;
  721. t[2] = s*x - 1.0;
  722. for( i=3; i<n; i++ )
  723.     t[i] = s * t[i-1] - t[i-2];
  724.  
  725. s = 0.0;
  726. for( i=0; i<n; i++ )
  727.     s += t[i] * ch[n-i-1];
  728. return(s);
  729. }
  730.